function [params_hat, covMat, params_ses, modelFits] = estimate_model(modelId, data, varargin)
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	% This function estimates the model by maximum likelihood (using Newton-Raphson algorithm).
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	%%%%% Inputs:
	% modelId;					integer
	% data:			           object:
	% 	.dims:				     	object:
	%		.dims1:						cell(NumParts,1)
	%			{ii}:						integer: gives dim1_ii that corresponds to Xparts{ii}
	%		.Xpart_2_NumX:				1 x NumParts: gives integers (dim2 of Xparts{ii}.X)
	%		.Xpart_2_NumX_FEs:			1 x NumParts: gives integers (dim2 of Xparts{ii}.X_FEs)
	%		.Xpart_2_Num_FE_vals:		cell(1,NumParts)
	%			{ii}:						1 x NumX_FEs_i  --> gives integer: number of possible values for Xparts{ii}.X_FEs(:,ff)
	%		.NumFEvals2Keep:			cell(1,NumParts) --> gives integer
	%		.NumParts
	%		.NumObs
	%		.NumParams
	%    .LL_offset:               1 x 1
	%	 .Xparts:				   cell(NumParts,1)
	%	 	{ii}:				       object
	%			.X:					       dim1_ii x NumX
	%			.X_FEs:					   dim1_ii x NumX_FEs
	%			.NumX_FE_vals:			   1 x NumX_FEs: gives integer
	%    .N:                       NumObs x 1    (if modelId==1)
	%    .logLambdaOffset:         NumObs x 1     (if modelId==2)
	%    .Y:                       NumObs x 1
	% varargin{1}=params0:		   NumParams x 1
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	%%%%% Outputs:
	% params_hat:				   NumParams x 1
	% covMat:					   NumParams x NumParams
	% params_ses:				   NumParams x 1
	% modelFits:				   object:
	%	.LL:							1 x 1
	%	.AIC:							1 x 1
	%	.BIC:							1 x 1
	%	.MAE:							1 x 1
	%	.MSE:							1 x 1
	%	.RMSE:							1 x 1
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	
	% Read dimensions
	NumParams = data.dims.NumParams;
	NumObs    = size(data.Y,1);
	disp(sprintf('Number of observations: %d', NumObs));
	disp(sprintf('Number of parameters: %d',   NumParams));
	
	% Read params or set default value
	if nargin >= 3
		params0 = varargin{1};
	else
		params0 = zeros(NumParams,1);
	end
		
	%%% Define objective function
	myfunc = @(params_val) loglikelihood(modelId, data, params_val, true);
	
	%%% Estimate by Newton-Raphson
	tic
	params_hat = Newton_Raphson(myfunc, params0, 1e-8, 1000, 1, true);
	toc
	
	%%% Calculate standard errors
	[LL,~,FisherInfo] = myfunc(params_hat);
	LL = -LL;
	covMat = inv(FisherInfo);
	params_ses = sqrt(diag(covMat)); % Standard errors
	
	%%% Compute AIC and BIC
	modelFits.LL = LL;
	modelFits.AIC = 2*NumParams - 2*LL;
	modelFits.BIC = NumParams*log(NumObs) - 2*LL;
	
	%%% Compute (static) prediction errors
	% Compute predictions
	dims = data.dims;
	params_parts = split_combine_parts(dims, 1, params_hat);
	V = compute_Xbeta2(dims, data.Xparts, params_parts);	% NumObs x 1
	if modelId == 1  % Logistic-Binomial --> Ypred = data.N .* p
		Ypred  = data.N./(1+exp(-V)); % NumObs x 1
	end
	if modelId == 2  % Poisson --> Ypred = lambda
		Ypred  = exp(data.logLambdaOffset + V); % NumObs x 1
	end
	clear V;
	
	% Evaluate prediction errors (RMSE, MAE, etc.)
	pred_errors = data.Y - Ypred; % NumObs x 1
	clear Ypred;
	MAE  = mean(abs(pred_errors));
	MSE  = mean(pred_errors.^2);
	RMSE = sqrt(MSE);
	clear pred_errors;
		
	%%% Output results to files
	modelFits.MAE = MAE;
	modelFits.MSE = MSE;
	modelFits.RMSE = RMSE;
end
